m  FILE  COPt 


AFOiR.TR.  on  ..  „  (•> 

O  9  -  0  0  89^— 


Annual  Technical  Report  #1 
November  1, 1987  -  November  1, 1988 


ANALYSIS  OF  FLOW-, 

THERMAL-  AND  STRUCTURAL-INTERACTION 

OF 

HYPERSONIC  STRUCTURES 
SUBJECTED  TO  SEVERE  AERODYNAMIC 

HEATING 


Contract  No,  F49620-C-0001 


Air  Force  Office  of  Scientific  Research 
Boilding  410 

Bolling  AFB,  DC  20332-6448 


ELECTE 
FEB  1  6 1989 


I 


ptmimmox  -statement 

approved  tot  public  teleaa 
Distribution  UoTobiUid 


rtMENT  Cl 

c  rbieaaej  1 
limited  i 


TR-88-12 


/sm  r~*  *' 


me  cpwwr*rtUf<At  wcatAmo  comfany.  me 


THE  COMPUTATIONAL  MECHANICS  CO.,  INC. 
_  3701  North  Lamar,  Suite  201 

Austin,  Texas  78705 
(512)  4674)618 


89_ 2  15  16 


REPORT  DOCUMENTATION  PAGE 


form  Approved 
QMBNo.  0704-Q188 


la.  REPORT  SECURITY  CLASSIFICATION 
Unclassified 


2a.  SECURITY  CLASSIFICATION  AUTHORITY 


2b.  DECLASSIFICATION,/ DOWNGRADING  SCHEDULE 


1b.  RESTRICTIVE  MARKINGS 

None 


I  *  OKTRIftLLTinN /AVAILABILITY  OF  REPORT 

Ipumfl  for  public  release . 
distribution  unlimited 


4.  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

TR-88-12 

• 

6a.  NAME  OF  PERFORMING  ORGANIZATION 
Computalonal  Mechanics 
Company,  Inc. 

6b.  OFFICE  SYMBOL 
(If  applicable) 

6c.  ADDRESS  (Oty,  Stite,  end  HP  Code) 

3701  North  Lamar,  Suite  7.01 
Austin  TX  78705 

5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 

AFOSR.TR.  8  9-  0  089 


■  r  o  ri 


8a.  NAME  OF  FUNOING/ SPONSORING 
ORGANIZATION  USAF,  AFSC 
Air  Force  Office  of  Scientifi 


8  c.  AOORESS  (City,  stare,  and  Z/P  Code) 
Bolling  Air  Force  Base 
DC  20332  ft  Id’ 


7b.  ADDRESS  (City,  State,  and  HP  Code) 

Su.d\-e. 

be. 

9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER  jj 

F49620-88-C-0001 

_ J 

10.  SOURCE  OF  FUNDING  NUMBERS  I 

(i  < 


1 1 .  TITLE  (inciude  Security  Classification) 

“Analysis  of  Flow-,  Thermal-,  and  Structural-Interaction  of  hypersonic 
Structures  Subjected  to  Severe  Aerodynamic  Heating” 


12.  PERSONAL  AUTHOR^} 

J,  T,  Oden  /  E,  A.  Thornton 


13a.  type  OF  REPORT  13b  time  COVERED  jU.  DATE  OF  REPORT  (Year, Month.  Dey)  IIS.  PAGE  COUNT 

Annual  Technical  from  11-1-87tq  I 1— I— November  30,  1988  | 


16.  SUPPLEMENTARY  NOTATION 


cosai'i  cooes 


GROUP 


November  30,  1988 


19.  A8ST 


tinue  on  reverse  if  necessary  snd  identify  by  block  number) 


This  first  annual  report  presents  progress  in  the  modelling  of  hypersonic 
fluid-therm ai-stnjeturai  interaction.  In  this  phase  of  the  effort,  the  basic  mechanisms  of 
heat  transfer  and  fluid-structure  interaction  were  identified.  Mathematical  models  for  heat 
transfer,  structural  deformation  and  fluid  flow  analysis  were  formulated.  New  unified 
viscoplastic  theories  were  (adapted  for  the  modelling  of  complex  visco-el asto-plastic 
structural  deformation,  with  temperature-dependent  material  properties.  A  general 
procedure  for  the  analysis/ of  fl uid •  thermal-structure  interaction  was  formulated,  and 
relevant  finite  element  cod^s  were  develop.  These  problems  were  applied  in  the  solution 
of  representative  examples  of  thermo-structural  analysis  of  structures  subject  to 
aerodynamic  heating.  <JInese  programs  of  adaptive  finite  clement  analysis  of  viscous  flow 
were  also  developed  and  validated  on  selected  examples.  Finally,  directions  for  further 
research  in  the  project  were  developed. 


20.  DISTRIBUTION /AVAILABILITY  OFACSTRACT 
6^ NC LASS IFIEO /UNLIMITED  □  SAME  AS  RPT 


224  NAME  OF  RESPONSIBLE  INDIVIDUAL 

Dr.  Anthony  K.  Amos 


DO  Form  1473.  JUN  86 


21.  ABSTRACT  SECURITY  CLASSIFICATION 
GHStic  USERS  Unclassified 


22b  telephone  (include  Are*  Code) 
(202)  767-4937 


Previous  editions  tre  obsolete. 


iBammmjma 

mL"  ....  .  .  IQ 


SECURITY  CLASSIFICATION  OF  THIS  PAG 


UNCLASSIFIED 


Annual  Technical  Report  #1 
November  1, 1987  -  November  1, 1988 


ANALYSIS  OF  FLOW-, 

THERMAL-  AND  STRUCTURAL-INTERACTION 

OF 

HYPERSONIC  STRUCTURES 
SUBJECTED  TO  SEVERE  AERODYNAMIC 

HEATING 


Contract  No.  F49620-C-0001 


Air  Force  Office  of  Scientific  Research 
Building  410 

Bolling  AFB,  DC  20332-6448 


jjAcCt'xieit  fi)t 

NTIS  CHAAI 

or  to  ua 

JiKjilotum _ 

o 

□ 

By 

Dst  it)  i|i>*  tj 


A 


Dist 


Codes 

Awl  ,jn d/or 


■me  COMfH/TATtO  HAL  lOCUAhUCS  COMPANY.  LNC. 


TR-88-12 


THE  COMPUTATIONAL  MECHANICS  CO.,  INC. 

3701  North  Lamar,  Suite  201 
Austin,  Texas  7870S 
(512)  467-0618 


Contents 


1  Introduction  1 

2  Research  Progress  3 

3  Thermal  Analysis  of  Hypersonic  Structures  6 

3.1  Engineering  Model  for  Coolant  Passages . 6 

3.2  Finite  Element  Formulation .  8 

4  Thermo— Viscoplastic  Structural  Analysis  10 

4.1  Initial  Value  Viscoplasticity  Problem .  11 

4.2  Constitutive  Model .  12 

4.3  Finite  Element  Formulation .  15 

4.4  Viscoplastic  Solution  Method .  16 

4.5  Variable  Time  Step  Algorithm  .  17 

5  Viscous  Flow  Analysis  19 

5.1  Problem  Formulation .  19 

5.2  Taylor- Galerkin  Algorithm .  21 

5.3  Aerothermal  Loads .  22 

6  Recent  Numerical  Results  24 

6.1  Thermo- Viscoplastic  Structural  Computations  . .  24 

6.1.1  Simplified  Model . 24 

6.1.2  Convectively  Cooled  Structure .  25 

6.2  Flow  Computations . 26 

6.2.1  Flow  Over  Flat  Plate  .  26 

6.2.2  Boundary  Layer  Flow . 27 

7  Research  Forecast  28 

8  References  29 


1  Introduction 


The  commitment  to  develop  the  National  Aerospace  Plane  (NASP)  has  generated  resur¬ 
gent  interest  in  the  technology  required  to  develop  structures  for  hypersonic  flight.  Such 
structures  will  be  exposed  to  aerodynamic  heating  of  unprecedented  magnitudes.  NASA 
estimated  surface  temperatures  for  NASP  are  shown  in  Fig.  1.  As  the  vehicle  acceler¬ 
ates  or  decelerates  at  hypersonic  speeds  in  the  earth’s  atmosphere,  shock  waves  will  sweep 
across  the  vehicle  and  interact  with  local  shocks  and  boundary  layers.  These  interactions 
introduce  severe  local  pressures  and  heating  rates.  A  recent  experimental  study  (ref.l) 
of  interacting  shock  waves  on  a  cylindrical  leading  edge  shows  heating  rates  ten  times 
undisturbed  levels. 

Leading  edges  of  engine  structures  present  a  signficant  design  problem  because  of  in¬ 
tense  local  heating  and  pressures.  Analysis  of  the  flow,  thermal  and  structural  behavior 
present  serious  computational  challenges  to  analysts  because  of  the  inherent  nonlinearities 
in  all  aspects  of  the  multi-disciplinary  problems.  Some  of  the  critical  computational  issues 
me  identified  in  ref.  2.  Critical  issues  include  the  difficulties  involved  in:  (1)  analysing  the 
viscous,  compressible  flow  and  predicting  the  lugh  local  aerodynamic  heating,  (2)  model¬ 
ing  and  analyzing  multi-mode  unsteady  heat  transfer  in  a  high  temperature  convectively- 
coolod  structure,  and  (3)  simulating  the  transient,  nonlinear  thermal-structural  response 
for  rapid  temperature  changes.  Preliminary  structural  analysis  of  an  impingement  cooled 
leading  edge  (ref. 3)  showed  high  local  plasticity  that  seriously  degraded  the  structure’s 
load  carrying  capacity  at  elevated  temperatures.  A  recent  therniostructural  analysis  with 
experimental  verification  (ref.4)  of  cowl  lip  designs  confirmed  that  significant  inelastic  ef¬ 
fects  occur,  In  the  experimental  study,  two  specimens  failed  due  to  burn-through  because 
of  intense  local  heating  or  because  of  loss  of  cooling. 

This  research  contract  addresses  fundamental  issues  relevant  to  the  development  of 
hypersonic  structures  subject  to  severe  thermal  loads.  The  specific  objectives  of  the  current 
research  include: 

1.  An  investigation  of  the  processes  and  theories  for  energy  transfer  between  lugh  tem¬ 
per  fhrrc  steady  hypersonic  flows  and  deformable  aerodynamic  structures. 
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2.  Identification  of  mechanisms  for  heat  transfer  in  hypersonic  structures  subjected  to 
intense  heating  rates.  Development  of  mathematical  models  for  high  temperature, 
inelastic  response  of  structures  considering  the  effects  of  time  dependent,  intense 
local  heating  and  pressures. 

3.  Development  of  a  unified  finite  element  mesh  refinement  scheme  for  interdisciplinary 
problems  including  flow,  thermal  and  structural  response. 

4.  Studies  of  unsteady,  transient  problems  of  hypersonic  flight  involving  structural  in¬ 
teraction  with  the  aerothermodynamics  of  the  external  flow. 


2  Research  Progress 


The  fluid- thermal-structural  problems  encountered  on  the  NASP  scramjet  engine  structure 
shown  in  Figure  2  are  illustrative  of  the  hypersonic  interaction  problems  being  considered. 
In  the  first  year  of  this  contract,  work  has  focused  on  the  physical  phenomena,  mathe¬ 
matical  formulations,  and  computational  methods  needed  to  solve  general  problems  of  this 
type. 

The  flow-thermal-structure  (F-T-S)  interaction  problem  being  considered  is  shown  in 
Figure  3.  Figure  3a  shows  a  low  temperature,  undeformed  elastic  structure  subjected  to 
low-level,  uniform  aerodynamic  heating.  In  this  state,  there  is  little  F-T-S  interaction.  In 
contrast,  Figure  3b  shows  a  structure  that  is  subjected  to  intense  localized  heating  char¬ 
acteristic  of  shock  interactions.  There  is  a  local,  growing  region  of  very  high  temperature 
plastically  deforming  material.  As  the  plastic  region  grows,  the  surface  deforms  into  the 
flow  producing  flow  recirculation.  The  flow  recirculation  intensifies  local  heating,  skin  fric¬ 
tion  and  surface  pressure.  As  the  heating  increases,  the  temperature  exceeds  the  melting 
temperature.  Thereafter,  a  more  complex  interaction  between  the  molten  material  and  the 
external  flow  occurs.  The  F-T-S  interaction  including  a  molten  region  is  not  understood; 
the  formulation  and  solution  of  this  problem  is  a  major  goal  of  the  research. 

As  steps  toward  achieving  tliis  goal,  the  following  objectives  have  been  accomplished: 

1.  A  formulation  of  a  coupled  F-T-S  interaction  problem  has  been  established  assuming: 

(a)  the  unsteady  heat  transfer  problem  is  described  by  conservation  of  energy  ne¬ 
glecting  mechanical  coupling  terms.  This  assumption  is  known  to  be  reason¬ 
able  for  the  materials,  deformation  rates  and  applied  heat  loads  encountered  in 
hypersonic  structures,  Nonlinear  heat  transfer  due  to  temperature  dependent 
material  properties  such  as  specific  heat  and  thermal  conductivity  arc  included 
as  well  as  surface  radiation 

(b)  the  structural  problem  is  assumed  to  be  quasi-static  and  viscoplastic.  The 
quasi-static  assumption  means  that  inertia  terms  in  the  momentum  equations 
may  be  neglected,  and  the  structural  problem  can  be  formulated  as  a  series 
of  equilibrium  problems  due  to  time-varying  temperature  computed  from  the 
unsteady  energy  equation.  The  vjseoplastic  model  assumes  that  strain  rates 
arc  separated  into  clastic  and  plastic  components.  The  plastic  strain  rate  com¬ 
ponent  is  represented  by  constitutive  equations  known  as  unified  viscoplastie 
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theories.  Several  unified  theories  exist  with  the  common  feature  that  strain 
rate  dependent  effects  such  as  plastic  flow,  creep  and  stress  relaxation  are  in¬ 
trinsic  to  the  formulation.  The  Bodner-Partom  theory  has  been  selected  because 
of: 

i.  the  availability  of  experimentally  determined  parameters  for  high  temper¬ 
ature  alloys,  and 

ii.  the  theory  has  been  validated  by  experiments  at  elevated  temperatures. 
The  viscoplastic  theory  assumes  that  temperatures  are  less  than  75  percent 
of  the  melting  temperature. 

(c)  The  flow  problem  is  governed  by  the  Navier-Stokes  equations  describing  com¬ 
pressible  flow.  The  principal  assumption  made  is  that  the  fluid  is  a  calorically 
perfect  gas.  This  assumption  is  legitimate  for  supersonic-hypersonic  flows  where 
the  Mach  number  is  less  than,  say,  seven.  Calorically  perfect  means  that  the 
gas  specific  heat  is  constant.  For  higher  Mach  numbers  where  the  gas  temper¬ 
ature  is  25G0K  or  higher  the  specific  heat  becomes  temperature  and  pressure 
dependent,  and  the  gas  becomes  chemically  reactive. 

2.  Three  finite  element  computer  programs  applicable  to  F-T-S  interactions  have  been 
placed  into  operation. 

(a)  A  family  of  finite  element  heat  transfer  programs  capable  of  computing  tem¬ 
peratures  in  Itigh  temperature  structures  is  now  operational.  These  programs 
model  nonlinear  conduction,  radiation  mid  forced  convection  heat  transfer  in 
conveetively  cooled  hypersonic  structures  at  elevated  temperatures. 

(b)  a  two'*diir.5msioasd  finite  element  thermo-viscoplastie  program  has  been  devel- 

oped  mn  is  based  on  the  Bodner-Partom  unified  viscoplastic  consti¬ 

tutive  model.  Highly  complex  rate-dependent  plasticity  and  creep  phenomena 
at  elevated  temperatures  can  be  simulated  by  the  program. 

The  use  of  the  Bodner-Partom  theory  in  a  finite  element  program  lias  significant 
advantages,  but  the  computational  effort  required  to  solve  for  tlie  state  variables 
in  the  constitutive  equations  is  extensive.  Sets  of  “stiff**  ordinary  differential 
equations  must  be  integrated  in  time.  Typically*  these  equations  require  very 
small  steps,  and  hence  large  computer  times  are  needed  to  capture  tlie  time- 
dependent  phenomena.  Time  steps  for  the  structural  respouse  are  typically 
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about  one-tenth  of  time  steps  required  for  the  transient  thermal  analysis.  This 
difference  in  time  steps  requires  special  techniques  for  interfacing  the  thermal 
and  structural  codes.  An  approach  to  interface  the  thermal  program  with  the 
viscoplastic  program  has  been  developed  and  implemented.  An  adaptive  time- 
step  algorithm  has  been  developed  and  implemented  to  give  reliable,  efficient 
solution  of  the  stiff  ordinary  differential  equations. 

(c)  A  two-dimensional  finite  element  viscous  flow  code  has  become  operational. 
The  program  solves  the  conservation  of  mass,  momentum  and  energy  equations 
for  compressible  viscous  flows.  The  conservation  equations  are  solved  using 
an  explicit  time-marching  algorithm.  An  adaptive  h  refinement /derefinement 
scheme  has  been  implemented.  A  “consistent”  method  for  computing  aerody¬ 
namic  surface  heating  rates,  skin  friction  (shearing  stresses)  and  pressures  has 
been  developed  and  implemented.  The  capability  of  the  viscous  flow  code  to 
predict  surface  quantities  is  under  evaluation. 

3.  Numerical  computations  have  begun  on  a  series  of  representative  problems,  Results 
of  these  computations  are  presented  later  in  this  report.  The  objectives  of  these 
computations  are  to; 

(a)  validate  the  codes,  aud 

(b)  study  the  basic  processes  that  characterize  F-T-5  interactions. 

4.  A  comprehensive  literature  search  focused  on  the  behavior  and  matlnnatical  models 
of  metallic  materials  at  elevated  temperatures  has  been  accomplished.  A  significant 
number  of  relevant  papers  and  reports  have  been  acquired. 

5.  Visits  to  Southwest  Research  Institute  and  Texas  A  k  M  University  have  been  made. 
Information  on  recent  research  in  unified  viscoplastic  mathematical  models  was  ob¬ 
tained  as  well  as  data  describing  material  behavior  at  elevated  temperatures.  Imple¬ 
mentation  of  these  theories  into  finite  element  codes  was  discussed. 


3  Thermal  Analysis  of  Hypersonic  Structures 


Convectively  cooled  structures  are  strong  candidates  for  use  in  hypersonic  flight  vehicles. 
For  hypersonic  flight,  some  leading  edges  and  panels  require  active  cooling  systems  to 
keep  structural  temperatures  within  acceptable  ranges.  The  internal  flow  in  the  coolant 
passage  has  a  predominant  role  in  the  thermal  response  of  a  hypersonic  structure  subject 
to  external  heating.  A  cross-section  of  a  typical  convectively  cooled  structure  is  shown 
in  Figure  4.  An  aerodynamic  skin  and  a  coolant  passage  with  internal  heat  exchanger 
protect  the  primary  structure  from  the  aerodynamic  heating.  The  thin,  typically  metallic, 
aerodynamic  skin  transfers  the  energy  of  the  aerodynamic  heating  to  a  low  temperature 
coolant  flow  through  the  heat  exchanger  fins  that  connect  the  aerodynamic  skin  to  the 
primary  structure.  In  a  typical  engine  structure,  the  coolant  is  cold  hydrogen  that  later  is 
used  as  the  propulsion  system  fuel. 

Heat  transfer  in  the  aerodynamic  skin  consists  of  conduction  combined  with  surface 
radiation.  Heat  transfer  between  the  aerodynamic  skin,  the  heat  exchanger  surfaces  and 
the  primary  structure  is  by  conduction  at  the  solid-fluid  interface.  The  dominant  mode  of 
heat  transfer  in  the  coolant  flow  is  forced  convection. 

The  representation  of  the  heat  transfer  in  the  coolant  passage  is  the  critical  step  in  the 
heat  transfer  analysis.  There  are  two  basic  representations  that  can  be  employed.  The 
first,  denoted  here  as  the  engineering  model,  is  based  upon  a  number  of  assumptions  that 
greatly  simplify  the  problem  into  a  single  energy  equation  with  a  specified  mass  flow  rate. 
Detailed  computation  of  the  fluid  velocity  components  and  temperatures  is  not  required. 
The  second  representation  idealises  the  coolant  flow  as  a  continuum  model,  and  the  par  tial 
differential  equations  describing  conservation  of  mass,  momentum  and  energy  are  solved 
simultaneously  to  obtain  fluid  velocity  and  temperature  distributions.  This  latter  model  is 
the  most  accurate,  but  it  is  also  considerably  more  expensive  than  the  engineering  model. 
As  a  first  step  towards  understanding  fiowdhei  mahstructural  interactions  characteristic  of 
hypersonic  structures  the  engineering  heat  transfer  models  employed. 

3.1  Engineering  Model  for  Coolant  Passages 

The  basic  features  of  the  engineering  model  can  be  developed  from  the  idealisation  shown 
in  Fig.  5.  A  segment  of  the  coolant  passage  of  width  u>  is  shown;  for  simplicity,  only  the 
upper  one-half  of  the  coolant  passage  with  the  aerodynamic  skin  is  shown.  The  engineering 
formulation  (ref.fi)  is  based  on  the  following  assumptions: 
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1.  The  thermal  energy  state  of  the  fluid  is  characterized  by  the  fluid  bulk  temperature 
Tp  which  varies  only  in  the  flow  direction,  i.e.  Tp(x,  t). 

2.  The  flow  is  represented  by  the  mass  flow  rate  rh  in  the  coolant  passage  specified  by 
rh  =  ppApVp  where  pp  is  the  coolant  density,  Ap  is  the  cross-sectional  area  of  the 
coolant  passage,  and  Vp  is  the  coolant  mean  flow  velocity. 

3.  A  convection  coefficient  h  is  defined  such  that  the  heat  flux  qp  transferred  between 
the  structure  and  the  coolant  may  be  expressed  as 

4  =  KTs  -  Tf) 


where  Ts(x,t )  denotes  the  structural  temperature  at  the  fluid-solid  interface. 


4.  The  convection  coefficient  h  may  be  expressed  as  a  function  of  the  fluid  bulk  tem¬ 
perature  alone  by  using  analytieal/empirical  equations  for  the  Nusselt  number, 

..  hD 

jVt*  K2  ~— 

kp 

where  D  is  the  hydraulic  diameter  of  the  coolant  passage,  and  kp  is  the  thermal 
conductivity  of  the  fluid  coolant.  The  convection  coefficient  as  well  as  the  solid  and 
fluid  thermal  parameters  may,  in  general,  be  temperature  dependent. 


With  these  assumptions,  energy  balances  on  the  aerodynamic  skin  and  coolant  give  the 
governing  conser  vation  equations. 


Fluid: 


Bn  ,  BTp,  .  OTp 


wh(Ts  ***  Tp)  4*  ppcpAp 


dt 


0 


4*  wh(T$  —  Tp)  +  4  p$GsA$~*j~  s*  4 


dx 


where  the  subscripts  F  and  5'  denote  the  fluid  and  solid,  resjjectively.  In  these  equations, 
c  denotes  specific  heat,  <f  is  the  Stefan- Boltzmann  constant,  and  (  is  the  surface  enussivUy. 
Because  of  the  temperature  dependence  of  the  thermal  parameters  mid  the  radiation  term, 
Fqs.  (3.1)  and  (3.2)  constitute  a  nonlinear  set  of  partial  differential  equations. 

The  heat  exchanger  fins  are  not  included  explicitly  in  this  model.  However,  the  heat 
transfer  between  the  heat  exchanger  fins  can  approximately  be  taken  into  consideration 


though  the  use  of  art  effective  width  w  which  represents  the  area  over  which  the  convective 
heat  exchange  occurs. 
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3.2  Finite  Element  Formulation 

A  typical  finite  element  representing  Eqs.  (3.1)  and  (3.2)  is  characterized  by  fluid  and 
fluid-solid  interface  nodes.  The  element  shown  in  Fig.  6  has  two  fluid  nodes  (I  and  J), 
and  two  fluid-solid  interface  nodes.  Within  the  element,  the  fluid  and  solid  temperatures 
are  express  as 

Tf  =  [N(x)]{TF} 

TS  =  {N(x)){Ts}  (3.3) 

where  (Ar(a:)]  are  the  element  interpolation  functions.  Following  usual  finite  element  pro' 
cedures  (ref.  o),  the  discritized  equations  for  an  element  of  length  L  may  be  derived  in  the 
form 


‘  CF  0 

*r  1 

< 

0  Cs 

l  Ts 

'*  +  Kk  +  A>  -A* 

■  i 

Tf  1 

-  k\  KS  4-  k‘H 

l  J*  . 

where  the  element  capacitance  matrices  are  given  by 
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and  the  element  conductance  matrices  are 
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The  element  equations  given  in  Eqs.  (3.4)  show  that  the  coolant  passage  model  can  be 
regarded  as  an  assembly  of  elements  where  each  element  represents  a  single  heat  transfer 
mode.  Thus  we  can  represent  the  coolant  passage  in  terms  of  two  convective  elements:  (1) 
a  mass-transport  element  (Fig.  7a),  and  (2)  a  surface  convection  element  (Fig.  7b).  The 
mass  transport  element  represents  the  downstream  convective  heat  transfer  due  to  mass 
how;  it  has  a  capacitance  matrix  [Cp]  as  well  as  conductance  matrices  [/<"„]  and  [ Kp ]•  The 
surface  convection  element  represents  the  convection  heat  exchange  between  the  coolant 
fluid  and  the  solid;  it  does  not  have  a  capacitance  matrix,  but  it  has  a  conductance  matrix 
[A’a]  linking  the  solid  and  fluid  nodes.  These  coolant  passage  elements  are  assembled  into 
the  finite  element  thermal  model  for  the  complete  convectively  cooled  structure  that  has 
conduction  elements,  radiation  elements,  and  surface  elements  for  convection  to  a  specified 
convective  exchange  temperature. 

The  unsteady  thermal  analysis  is  nonlinear  because  of  temperature  dependent  thermal 
properties  and  surface  radiation.  The  equations  are  solved  by  time  marching  with  the 
Crank-Nicolson  algorithm;  at  each  time  step,  the  nonlinear  algebraic  equations  are  solved 
by  Newton-Raphson  iteration. 

.  _  > 

i  -  ' 
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4  Thermo-Viscoplastic  Structural  Analysis 

Unified  viscoplastic  constitutive  models  have  evolved,  over  the  last  twenty  years,  to  provide 
means  for  analytically  representing  material  response  from  the  elastic  through  the  plastic 
range  including  strain-rate  dependent  plastic  flow,  creep  and  stress  relaxation.  The  theo¬ 
ries  are  guided  by  physical  considerations  including  dislocation  dynamics  and  are  based  on 
the  principles  of  continuum  mechanics.  The  first  multi-dimensional  formulation  of  elastic- 
viscoplastic  constitutive  equations  was  due  to  Bodner  and  Partom.  Since  then  a  number 
of  constitutive  models  have  appeared;  many  of  these  theories  are  summarized  in  review 
articles  that  appear  in  reference  6.  A  NASA-Lewis  sponsored  research  program  (HOST) 
conducted  by  the  Southwest  Research  Institute  recently  concluded  a  four  year  research  ef¬ 
fort  (ref.  7-8)  to  further  develop  unified  constitutive  models  for  isotropic  materials  and  to 
demonstrate  their  usefulness  for  analysis  of  high  temperature  gas  turbine  engines.  One  re¬ 
sult  of  this  study  is  material  property  data  for  nickel-based  alloys  over  a  wide  temperature 
range.  The  unified  models  employed  were  those  of  Bodner-Partom  and  Walker. 

Unified  viscoplastic  theories  have  been  implemented  by  a  number  of  finite  element 
researchers.  Under  the  NASA  HOST  program,  the  Walker  model  was  implemented  in 
the  MARC  finite  element  program  (ref.  9)  and  used  to  analyze  the  thermo-viscoplastic 
response  of  a  turbine  blade  under  simulated  flight  conditions.  In  another  recent  finite 
element  application  (ref.  10),  the  Bodner-Partom  and  Walker  theories  were  compared  for 
a  thin  circular*  plate  subject  to  highly  localized,  transient  heating. 

In  this  research  the  Bodner-Partom  constitutive  model  is  employed,  and  the  finite  el¬ 
ement  approach  developed  in  reference  11  for  the  isothermal  case  is  extended  to  include 
thermal  effects.  The  behavior  of  a  thermo-viscoplastic  structure  subjected  to  aerodynamic 
heating  is  analyzed  assuming  that:  (1)  thermo-mechanical  coupling  in  the  conservation 
of  energy  equation  can  be  neglected,  (2)  the  structural  response  is  quasi-static,  and  (3) 
deformations  are  infinitesimal.  With  these  assumptions,  an  unsteady  thermal  analysis  may 
be  performed  fust  to  determine  the  temperatures.  Then,  using  these  temperatures,  the 
structure’s  viscoplastic  response  is  determined.  The  solution  is  thus  obtained  by  sepa¬ 
rately  solving  initial  boundary  value  problems  for  first  the  thermal  and  then  the  structural 
response. 
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4.1  Initial  Value  Viscoplasticity  Problem 

Consider  a  viscoplastic  structure  occupying  a  region  Q.  with  boundary  The  behavior 
of  the  structure  is  described  by  the  following  system  of  differential  equations: 

1.  Equilibrium  in  rate  form, 

crijj  +  i>i  =  o  (4.1) 

where  cr;j  is  the  stress  tensor,  are  the  body  forces  per  unit  volume,  and  the  sum¬ 
mation  convention  is  employed. 

2.  Kinematic  relation  for  velocity  gradients, 


Qj  -  «§  +  *fj  —  +  Uj.i) 


(4.2) 


where  e,j  is  the  total  strain  and  superscripts  E  and  P  denote  elastic  and  inelastic 
strain  components,  respectively.  The  displacement  rates  are  ttj. 

3.  Constitutive  relations 

&{j  ss  Eijkiiitf  —  E\jki o fci AT 

2k)  (no  sum)  (4.3) 

-  <Ji{<rij,Zk)  (no  sum) 

where  E^ta  is  Hooke's  tensor  of  elasticity  parameters,  and  aty  is  a  tensor  of  thermal 
expansion  parameters,  and  AT  represents  the  rate  of  the  change  in  temperature  from 
a  reference  temperature.  Both  Eiju  arid  aw  are  temperature  dependent.  The  con¬ 
stitutive  functions  are  /y  and  gt  where  Zt,  represents  internal  state  variables.  These 
functions  and  state  variables  characterize  the  viscoplastic  response  of  the  material. 


The  description  of  the  problem  is  completed  by  prescribing  the  boundary  and  initial 
conditions, 

ti  s  5 i  on 

(4.4) 

ctijtij  £=  S'j  on  d$h 
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where  a,-  are  prescribed  surface  displacement  rates,  rij  are  the  components  of  a  unit  nor¬ 
mal  vector,  and  a,-  are  prescribed  surface  traction  rates.  The  initial  conditions  include 
specifying  the  displacements,  stresses  and  internal  state  variables,  i.e. 

u,(x,  0),ati(x,  0),2,(x,0). 


4.2  Constitutive  Model 


We  now  describe  the  explicit  forms  (ref.  6-7)  of  the  Bodner-Partom  constitutive  equations 
used  in  the  thermo-viscoplastic  analysis.  The  Bodner-Partom  constitutive  model  is  of  the 
internal  state  variable  type  that  is  a  phenomenological  theory  based  on  the  physical  con¬ 
cepts  related  to  dislocation  dynamics.  The  model  has  gone  through  several  modifications 
and  extended  for  anisotropic  work  hardening  materials.  The  current  model  also  includes 
temper aure  effects. 


1.  Flow  Law, 

For  the  inelastic  strain  rate  component,  the  isotropic  form  of  the  Prandtl-Reuss  law 
is  assumed 


<5  =  «« 

«&  -  0  A>0 


(4.5) 


where  are  the  deviotoric  stress  components  given  by 


“  ^§ij*kk 


(4.6) 


5=s  0  denotes  plastic  incompressibility.  _ 

2.  Kinematic  Equations, 

Squaring  Eq.(4.5)  leads  to 

A2  =  Df/Jt  (4.7) 

where  D%  and  J-.  are  the  strain  rate  and  deviatoric  stress  invariants.  The  deviatoric 
stress  invariant  is 

(4.S) 
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The  relation  governing  inelastic  deformations  is  the  “kinetic  equation”,  and  the  form 
taken  by  Bodner-Partom  is 

DZ  =  D20exp[-(Z2/ZJ2)']  (4.9) 

where  D0  is  the  limiting  strain  rate  in  shear,  n  is  a  material  constant,  and  Z  i3 
interpreted  as  a  load  history  dependent  parameter  (herein  called  the  internal  state 
variable)  that  represents  the  hardened  state  of  the  material  with  respect  to  resistance 
to  plastic  flow.  Combining  Eqs.(4.5),  (4.7),  and  (4.8)  gives 

*5  =  [-|(2V3/1)”]  (4.10) 

3.  Evolution  Equations  of  Internal  State  Variable, 

The  internal  state  variable  Z  consists  of  isotropic  and  directional  components, 

Z~Zl  +  ZD  (4.11) 

The  evolution  equation  proposed  for  the  isotropic  hardening  component  (ref.  6)  is 

z'(i)  =  mxlz,  -  Zi(t)]Wp(t)  -  AiZ\  ZltthJk  (4.12) 
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with  the  initial  condition,  Z\ 0)  =  Za.  In  the  first  term,  Z\  is  the  limiting  (saturation) 
value  of  Z1 ,  m\  is  the  hardening  rate,  and  the  plastic  work  rate  is 


Wp  =  <j„t' 


(4.13) 


which  is  taken  as  the  meas\ire  of  hardening.  Z2  is  the  minimum  value  of  Z1  at  a 
given  temperature,  and  A\  and  r*  are  temperature  dependent  material  constants. 

The  evolution  form  of  the  directional  hardening  component  (ref.6)  is  defined  as 


Z°(f)  =  *«“«(<) 


(4.14) 


where  w,j  are  the  direction  cosines  of  the  current  stress  state, 


(4.15) 


The  evolution  equation  for  4(0  has  the  same  general  form  as  that  for  isotropic 
hardening  but  has  tensorial  character, 


4(0  «  mfauijit)  -  WWp(t) 


(4.16) 


where 

».;(<)  =  /MO/l/MOAKOI*  (4.17) 

and 

4(0)  —  0  (4.  IS) 

As  in  Eq.(4.12),  m2  is  the  hardening  rate.  „43  and  r2  are  temperature  dependent 
material  constants. 
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Other  variants  of  the  isotropic  and  directional  hardening  variables  are  described  in 
re£s.(6-8). 

4.3  Finite  Element  Formulation 

In  this  section,  the  finite  element  formulation  for  the  initial  value  thermo-viscoplasticity 
problem  (Section  4.1)  is  described.  The  development  follows  the  approach  of  Ref.  11,  and 
the  details  of  the  weak  formulation  are  presented  there.  Since  the  discussion  presented  in 
Ref.  11  is  for  the  isothermal  case,  this  formulation  extends  the  previous  approach  by  the 
inclusion  of  temperature  effects. 

The  finite  element  approach  approximates  displacement  rates  within  an  element  by 
taking  the  displacement  rates  as 

«  =  M{«}  (4.19) 

where  [JV]  are  the  interpolation  functions,  and  {6}  represent  the  nodal  displacement  rates. 
Using  the  strain-displacement  equations  in  rate  form,  Eq.  4.2,  an  element’s  strain  rates 
can  be  computed  from  Eq.  4.19  as 


{.)  =  [Z>l{i}  (4.20) 

where  (R]  is  the  strain-displacement  matrix.  Following  usual  finite  element  procedures, 
we  can  form  the  finite  element  equations  for  a  typical  element  as, 

(/C(T)]^} » {F,}  4-  {Fr}  4*  {&}  4-  {4}  (4.21) 


where  [K(T) )  is  the  element  stiffness  matrix,  and  the  terms  on  the  right-hand  side  are 
eleirent  load  vectors  due  to  the  rate  of  plastic  strains,  temperature,  surface  tractions  and 


body  forces,  respectively.  These  matrices  are  defined  by 


IWOI  =  JjBfwmda 

(4,22) 

{4)  = 

(4.23) 
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«■}  =  JalB\T[E(T)]{a(T)}ATdn 

(4.24) 

{&}  =  / 

Jd  n« 

(4.25) 

{ii}  =  f\N]T{b)da 

(4.26) 

where  Qe  denotes  the  element  volume,  and  dile  denotes  an  element  surface  where 
tractions  are  defined. 

The  temperature  affects  the  viscoplastic  structural  analysis  directly  in  three  ways:  (1) 
the  elasticity  matrix  (F(T)]  and  the  coefficients  of  thermal  expansion  {a(T)}  depend  on 
temperature,  (2)  nodal  loads  {Fr}  depend  on  the  local  temperature  rates,  and  (3)  several 
parameters  (see  section  4.2)  in  the  Bodner  Partom  constitutive  model  are  temperature 
dependent. 

4.4  Viscoplastic  Solution  Method 

Since  the  present  approach  uses  an  uncoupled  (sometimes  called  one-way  coupled)  for¬ 
mulation,  the  thermal  problem  is  solved  first  followed  by  the  viscoplastic  analysis.  The 
transient  thermal  problem  is  solved  by  time  marching  with  a  time  step  Air,  and  the  nodal 
temperatures  at  successive  times  ti,  #a,  *  *  •  are  obtained.  These  temperature  vectors  are 
used  as  input  to  the  structural  analysis. 

The  fust  computation  in  the  structural  analysis  is  to  solve  an  initial  statics  problem  if 
the  initial  temperature  distribution  T(x,0)  is  not  equal  to  a  uniform  reference  temperature. 
The  results  of  tins  analysis  are  the  initial  conditions  (displacements  and  stresses)  for  the 
transient  viscoplastic  analysis. . 

The  viscoplastic  analysis  time-marches  with  a  time  step  A f,.  Experience  has  shown 
the  time  step  required  for  the  structural  analysis  is  smaller  than  for  the  thermal  analysis, 
i.e.  At,  <  At}*.  At  intermediate  times  in  the  structural  analysis,  the  temperatures  are 
linearly  interpolated  from  the  temperatures  known  at  the  beginning  and  end  of  the  larger 
thermal  time  intervals. 

The  strategy  employed  in  the  viscoplastic  algorithm  is  as  follows:  with  the  initial  distri¬ 
bution  of  stress,  temperature  and  internal  var  iables  specified  use  the  equilibrium  condition 
(Eq.  4.21)  to  obtain  the  nodal  displacement  rates.  Then  integrate  the  constitutive  equa¬ 
tions  forward  in  time  at  the  element  Gauss  integration  points.  With  updated  values  of 


16 


the  stress,  temperature  and  internal  variables  at  the  new  time,  the  equilibrium  equation 
is  solved  again.  This  sequence  of  determining  the  nodal  displacement  rates,  then  advanc¬ 
ing  the  constitutive  equations  in  time  is  continued  until  the  desired  history  of  the  initial 
boundary-value  problem  has  been  obtained. 

Thus,  the  algorithm  proceeds  through  the  following  steps: 

1.  At  time  t,  initialize  for  each  element; 

2.  Calculate  Zk)  for  each  element; 

3.  Assemble  and  solve  [JC]{<5}  =  {F}  ; 

4.  Calculate  e,j  for  each  element,  {  e}  =  [iJ]{<$}; 

5.  Calculate  oq,  for  each  element,  {a}  =  [22]  {e  —  ep}  —  [2?]aAT; 

6.  Calculate  Z,  for  each  element,  Z,  =  Zk); 

7.  Integrate  <r<j,  Z  forward  for  each  element  to  get  <7,j  and  Z{  at  t  +  A t»\ 

8.  If  t 4-  At,  <  t fiUai  go  to  2,  otherwise  stop. 

The  computational  method  above  has  been  presented  for  a  constant  time  step  At,. 
Computational  experience  by  several  investigators  (see  Ref.  S,  10-11)  indicates  that  a  very 
small  time  step  can  be  required  because  of  the  “stiff”  nature  of  the  ordinary  differential 
equations  describing  the  internal  state  variables.  To  gain  improved  efficiency  and  reliability 
a  variable  time  step  algorithm  has  been  implemented  and  will  be  described. 

4.5  Variable  Time  Step  Algorithm 

In  the  constant  time  step  viscoplastic  algorithm  the  internal  state  variables  are  advanced 
in  time  with  the  conditionally  stable  Euler  forward  difference  algorithm.  The  variable  time 
step  algorithm  is  a  modified  Euler  scheme  using  a  truncation  error  criterion  (Ref.  12)  to 
adjust  the  time  step. 

For  simplicity,  consider  the  single  ordinary  differential  equation, 


V  ■-/(»,*)  (4.27) 
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The  solution  is  advanced  using  a  predictor-corrector  scheme.  The  predictor  phase  consists 
of  an  Euler  step: 

yt+At  =  yt  +  &tyt  (4.28) 

Vt+Ai  =  fiVt+At^  +  At)  (4.29) 


An  error  indicator  E  (Ref.  12)  is  then  computed  from 


1  At(yg.A<  -  y«)| 
2|y/+Atl 


(4.30) 


The  error  indicator  is  next  compared  with  a  preset  error  criterion,  and  if  the  criterion  is 
met,  the  time  step  is  sufficiently  small  enough  to  proceed  to  the  corrector  stage.  Otherwise, 
the  predictor  phase  Eqs.  (4.28-4.29)  is  repeated  with  a  smaller  time  step. 

The  corrector  phase  is  the  modified  Newton  scheme, 

V*vg  ’t’Vi+At)/ 2  (4.31) 

^  y*  Aij/avj  (4.32) 

A  flow  chart  depicting  the  adaptive  scheme  is  shown  in  Fig.  8.  The  flow  chart  shows 
how  the  time  step  is  either  reduced  or  increased  depending  on  the  error  indicator,  Eq. 
(4.30). 
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5  Viscous  Flow  Analysis 


Finite  element  analysis  of  compressible  flows  is  a  relatively  new  development,  but  significant 
progress  has  been  made  in  the  last  few  years.  Computational  techniques  were  developed 
first  for  inviscid  flows,  and  then  the  techniques  were  extended  for  viscous  flows.  The  first 
successful  algorithms  were  explicit  schemes;  current  research  is  focused  on  the  development 
of  implicit  schemes.  In  this  report  we  describe  an  explicit  finite  element  algorithm  (ref.  13) 
that  is  used  to  solve  the  Navier-Stokes  equations  for  the  external  flow.  A  basic  objective  of 
the  flow  analysis  is  to  predict  the  aerothermal  loads  acting  on  the  structure,  i.e.  the  surface 
heating  rate,  the  skin  function  (shearing  stress),  and  pressure  needed  for  the  thermal- 
structural  analysis. 

5.1  Problem  Formulation 

A  laminar  compressible  flow  with  a  calorically  perfect  gas  is  assumed.  The  Navier-Stokes 
equations  are  written  in  the  conservation  form, 


W)  .  d{E,-Ev)  d{F,-Fv) 

ir+ — & — + — §r —  (5,1) 

where  {17}  is  the  vector  of  conservation  variables;  {F/}  and  {F/}  are  inviscid  flux  compo¬ 
nents;  {Fv}  and  {Fv}  are  viscous  flux  components.  These  vectors  are  given  by 

{E/}r  =  [p  pu  pv  pEt 1 

{F/}r  -  [pu  pu2  +  P  puv  (pEt  +  P)u\ 

{F/}2’  s  [pv  pvu  pv 2  +  P  ( pEt  +  P)u]  (5.2) 

-  (0  rtu  uaxx  +  -  qx] 

{Fv}3,  -  (0  t**v  <rw  utzy  +  vcryv  -  gp] 

where  p  is  the  fluid  density,  u  and  v  arc  velocity  components,  Et  is  the  total  energy;  axx,  avv 
and  rXy  are  stress  components;  and  qxt  qv  are  heat  fluxes.  In  the  inviscid  flux  components 
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the  pressure  P  is  defined  for  a  perfect  gas  as 


P  =  (7  -l)p[E,  -(u’  +  v!)/2) 


(5.3) 


where  7  is  the  ratio  of  specific  heats.  In  the  viscous  terms,  the  stress  components  are 

ICldstCU  SO  VtiWvifc^  ~v 


<r„  =  !h(T)(2^-§) 

tw  =  p(D( $  +  &)  ■<«> 

»»*  =  §/‘(2’)(2|j  - 1|) 


where  T  is  temperature.  The  heat  fluxes  are  computed  by  Fourier’s  law 

*  =  *coSS 

«v  68  -&Cr)fy 


(5.5) 


The  temperature  dependent  viscosity  is  computed  from  Sutherland’s  law,  and  the  thermal 
conductivity  is  computed  assuming  a  Praadtl  number,  Pr  »  0.72. 

Equation  (5.1)  is  solved  subject  to  appropriate  initial  and  boundary  conditions.  Con¬ 
sider  a  domain  il  bounded  by  a  surface  da.  The  initial  conditions  consist  of  specifying 
the  distributions  of  the  conservation  variables  at  time  zero.  Typical  boundary  conditions 
may  include:  (1)  specifying  all  conservation  variables  for  supersonic  inflow  on  a  segment 
of  the  boundary  dSlu  (2)  requiring  that  $  •*  «  0  on  a  symmetry  plane,  B!h  (V  is  the 
velocity  vector,  and  n  is  a  unit  normal  vector)  and  specifying  the  heat  flux  normal  to  the 
surface  flflabe  zero,  (3)  using  appropriate  boundary  conditions  on  a  outflow  surface 
and  (4)  specifying  V  -  0  and  the  surface  temperature  T  on  aerodynamic  surfaces,  da4. 
On  supersonic  outflow  surfaces,  the  finite  element  formulation  supplies  appropriate  natural 
boundary  conditions  consisting  of  surface  integrals  of  flux  components. 


5.2  Taylor-Galerkin  Algorithm 

For  simplicity  the  solution  algorithm  will  be  given  for  the  single  scalar  equation, 


|+|ce ,-*)+£<*-*) -0 


(5.6) 


where  the  variables,  u,  Ei,  Ft,  Ev,  and  Fv  are  analogous  to  the  corresponding  vector  quan¬ 
tities  in  Eq.  (1).  Let  {u}n  denote  the  element  nodal  values  of  the  flow  variable  u(x,y,t)  at 
time  tn.  The  time  step  A t  spans  two  typical  times  tn  and  tn+l  in  the  transient  response. 
The  computation  proceeds  through  two  time  levels  <n+^  and  tn+\.  At  time  level  fn+|, 
values  for  u  that  are  constant  within  each  element  are  computed  explicitly.  At  time  level 
tn+i,  the  constant  element  values  are  used  to  compute  nodal  values  for  u. 


Time  Level  tn+j 

A  constant  element  value 


is  computed  from 


.4<*  ■«  /  [N)dA{u)»  -  f 


j  [f  ]  4ME,  -  EVr 


(5.T) 


where  A  denotes  the  area  of  element  and  (iVj  are  the  element  interpolation  functions. 
Quantities  such  as  {£/}*  represent  nodal  values  attime  t*. 

Time  Level  tH+.| 

Nodal  values  are  computed  from 


lA/llu)”*1  =  +  W**  +  +  W  +  W  (8.8) 


where 

\M\  =  Ja  WmdA  (5.3) 
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and 


{i?.2}n+*  =  {JV}d3[/(i;/)^+m(Fi),nH] 

w>n  =  -Ai  [L  {f }  w + L  { f}  w“  w 

w"  =  At  f  {N)iN]is[HEV}:+m{Fv}:\ 

Jou* 


(5.10) 

(5.11) 

(5.12) 

(5.13) 


In  Eqs.  (5.11)  and  (5,13),  /  and  m  are  the  components  of  a  unit  vector  normal  to  the 
boundary.  The  load  vectors  {Ri}  and  {11+}  represent  natural  conditions  and  are  evaluated 
on  outflow  surfaces. 

The  mass  matrix  (M),  Eq.  (5.9),  is  diagonalized  to  yield  a  “lumped  mass  matrix” 
producing  an  explicit  scheme  where  the  new  nodal  values  can  be  computed  directly  from 
Eq.  (5.8). 

5.3  Aerothermal  Loads 

Accurate  computation  of  surface  quantities  such  as  heat  duxes  is  an  important  objective  of 
the  viscous  flow  analysis.  By  Fourier’s  law,  heat  duxes  may  be  computed  from  temperature 
derivatives  at  the  fluid/solid  interface,  in  typical  elements*  temperatures  vary  linearly,  and 
temperature  gradients  normal  to  the  interlace  are  constant.  Heat  fluxes  computed  front 
these  derivatives  are  not  suflidently  accurate,  arm  conservation  is  not  guaranteed.  As 
an  alternate  approach,  a  finite  element  scheme  based  on  Eq.(S.S)  is  used.  Tins  latter 
approach  does  not  require  the  computation  of  temperature  derivatives,  and  conservation 
is  guaranteed. 

VVe  will  illustrate  the  approach  for  surface  heat  fluxes,  but  skin  friction  can  be  computed 
in  a  similar  fashion.  For  heat  fluxes,  wo  use  the  conservation  of  energy  equation  assuming 
the  computation  lias  reached  steady  state,  i.e.  {u}**1  —  {«}*.  On  no-slip  boundary 
surfaces  where  normal  heat  fluxes  are  to  be  computed,  Eq.  (5.8)  reduces  to 

ir}}  +  W^0 


since  the  inviscid  fluxes  in  the  energy  equation  are  zero  on  the  no-slip  boundary.  Substi¬ 
tuting  for  {R3}  using  Eq.  (5.10)  and  evaluating  {i?4}  at  the  surface  produces 


JjN}[N\da{tn] 


[. N]dA{Ev } 
[N}dA{Fv} 


(5.14) 


where  {?n}  are  nodal  heat  fluxes  for  an  element.  To  aid  in  the  computation,  the  coefficient 
matrix  for  {<?„}  in  this  equation  can  be  diagonalized  to  produce  an  explicit  scheme  for 
directly  computing  the  nodal  heat  fluxes. 


6  Recent  Numerical  Results 


Presented  here  are  some  recent  numerical  results  illustrating  progress  in  code  develop¬ 
ment  and  progress  in  understanding  basic  F-T-S  interaction  phenomena  for  hypersonic 
structures. 

6.1  Thermo-Viscoplastic  Structural  Computations 

6.1.1  Simplified  Modal 

To  gain  a  preliminary  understanding  of  the  behavior  of  a  convectively  cooled  structure 
subjected  to  convective  heating,  a  simple  ID  bar  model  was  investigated.  A  segment  (Fig. 
9a)  of  the  aerodynamic  skin  is  modeled  assuming  uniform  temperature.  The  coolant  is 
assumed  to  have  a  constant,  specified  temperature,  and  the  skin  is  heated  convectively. 
A  strong  time  valuation  of  the  convection  coefficient,  (Fig.  9b)  simulates  the  passage  of 
a  moving  shock.  The  resulting  transient  temperature  and  bar  compressive  stress  exhibit 
features  to  be  expected  in  a  more  complex  model. 

The  transient  temperature  history  (Fig.  10a)  shows  a  rapid  rise  following  the  sudden 
increase  in  the  convective  coefficient,  and  then  a  smooth  decay  after  the  large  convective 
heating  is  removed.  The  convective  cooling  causes  the  temperature  to  return  quickly  to 
equilibrium.  The  initial  rapid  increase  in  temperature  causes  the  bar  to  yield  in  com¬ 
pression  very  early  in  the  response  (Fig.  10b).  After  the  temperature  begins  to  decay 
at  t  —  0.5s,  the  stress  rapidly  changes  to  tension  hence  as  the  temperature  returns  to 
equilibrium  a  large  residual  tensile  stress  remains  in  the  bar.  The  tensile  stress  is  induced 
because  the  initial  high  temperature  causes  the  bar  to  expand  which  induces  compressive 
yielding.  Tliis  yielding  tends  to  permanently  shorten  the  bar.  However,  the  bar’s  fixed 
end  boundary  conditions  prohibit  the  bar  from  changing  in  length  so  that  a  tensile  stress 
develops  as  the  material  cools. 

In  the  initial  viscoplastic  analysis  a  fixed  time  step  of  0.001  s  was  used,  and  1200 
time  steps  were  required  for  the  analysis.  The  variable  time  step  algoritlnn  was  then 
implemented,  and  the  problem  was  resolved.  Fig.  10c  shows  the  history  of  the  variable 
time  step  and  indicates  that  the  new  analysis  required  only  213  steps  —  a  substantial 
savings.  The  figure  shows  that  in  the  MUat"  part  of  the  stress  response  a  large  time  step 
was  used,  but  near  t  —  0  and  again  at  t  =  0.5s  when  the  stress  is  changing  rapidly,  small 
time  steps  are  needed  to  capture  the  response  accurately.  This  example  shows  that  the 
adaptive  time  step  algorithm  is  an  important  step  towards  producing  an  efficient,  accurate 
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solution. 


6.1.2  Convectively  Cooled  Structure 

A  more  realistic  2D  model  of  a  convectively  cooled  hypersonic  structure  is  shown  in  Fig. 
11.  The  model  represents  a  segment  of  a  convectively  cooled  structure  such  as  a  wall 
of  a  scramjet  engine  fuel  injection  strut.  The  finite  element  thermal  model  (Fig.  11a) 
includes:  (1)  conduction  heat  transfer  in  the  aerodynamic  skin  and  primary  structure,  (2) 
convective  heat  transfer  between  the  walls  of  the  coolant  passage  and  coolant,  (3)  mass 
transport  convection  in  the  coolant  which  has  an  unknown  bulk  temperature,  and  (4) 
surface  radiation  on  the  aerodynamic  skin.  The  aerodynamic  skin  is  uniformly  convectively 
heated  over  its  length,  and  superimposed  on  the  uniform  heating  is  a  local,  intense  heating 
simulating  a  transient  shock.  The  transient  heating  is  induced  by  the  time-dependent 
convection  coefficient  shown  in  Fig.  lib.  Thermal  properties  and  data  are  tabulated  in 
Table  1.  The  properties  are  representative  of  a  high  temperature  nickel  superalloy. 

In  the  plane  strain  finite  element  model,  the  primary  structure  and  aerodynamic  skin 
have  unit  thickness,  but  the  heat  exchanger  in  the  coolant  passage  has  a  thickness  of  0.1 
in.  to  provide  an  effective  structural  stiffness.  The  wall  segment  has  fixed  displacement 
boundary  conditions  at  the  left  and  right  ends.  Data  for  the  Bodner-Partom  viscoplastic 
model  are  tabulated  in  Table  2. 

The  thermal  response  of  the  wall  segment  is  shown  in  Fig.  12.  Figure  12a  shows  the 
time  history  of  a  point  on  the  aerodynamic  skin  directly  under  the  transient  heating;  Fig. 
12b  shows  contours  of  temperatures  at  t  —  0.5s  when  temperatures  are  maximum.  The 
temperature  history  is  qualitatively  similar  to  the  results  obtained  in  the  ID  model  showing 
the  rapid  skin  temperature  rise  and  fall  with  the  simulated  shock  heating.  The  temperature 
contours  show  relatively  steep  thermal  gradients  induced  by  the  local  heating,  and  that  the 
coolant  substantially  limits  the  extent  of  the  induced  high  temperatures.  Thus  the  high 
temperatures  are  confined  to  the  aerodynamic  skin,  and  the  primary  structure  experiences 
only  small  temperature  changes.  The  temperature  gradients  in  the  skin  particularly  at 
the  coolant-skin  interface  are  not  predicted  with  high  accuracy  because  of  the  engineering 
model  of  the  coolant  heat  transfer.  However,  local  temperature  levels  in  the  skin  are 
reasonably  accurate  since  net  energy  transfer  to  the  coolant  is  modeled  satisfactorily. 

Stress  histories  at  three  points  through  the  skin  thickness  are  shown  in  Fig.  13.  The 
stress  histories  follow  the  temperature,  and  under  the  intense  local  heating  stresses  are 
very  similar  to  the  results  obtaiued  from  the  ID  model.  At  this  location  the  skin  yields 
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through  most  of  its  thickness.  After  the  heating  ceases  there  is  a  rapid  decay  of  stress,  and 
under  the  intense  local  heating  there  are  residual  tensile  stresses. 

The  viscoplastic  stress  histories  are  compared  with  stresses  predicted  assuming  elastic 
behavior  in  Fig.  14.  The  elasticity  assumption  is  clearly  not  permissable  because  predicted 
stresses  are  too  high,  and,  of  course,  no  residual  stresses  are  predicted.  Deformations  (not 
shown)  predicted  by  the  elastic  model  are  also  too  small. 

The  two  dimensional  character  of  the  stress  components  ax  and  cry  are  shown  in  Figs.  15 
and  16,  respectively.  The  stress  distributions  are  shown  at  three  times  (t  =  0.1s,  0.5s,  1.0s) 
in  the  response.  The  stresses  at  the  t  =  1.0s  are  the  residual  stresses.  The  figures  show 
that  in  addition  to  the  high  ax  stresses  induced  in  the  skin  that  cause  yielding,  there  are 
also  high  stresses  in  the  heat  exchanger  region.  Although  the  structural  model  of  the  heat 
exchanger  region  is  crude,  high  tensile  stresses  ay  can  potentially  cause  a  bond  failure  at 
the  heat  exchanger/skin  joint.  Such  a  failure  would  cause  a  loss  of  heat  transfer  between 
the  skin  and  heat  exchanger  resulting  in  excessively  high  skin  temperatures.  This  is  a 
possible  failure  mode  that  needs  further  investigation. 

The  deformed  structure  and  principal  plastic  strains  at  t  —  0.5s  are  shown  in  Fig. 
17  and  18,  respectively.  Peak  deformations  are  relatively  small,  but  the  effect  of  internal 
coolant  pressure  was  neglected.  The  effect  of  this  pressure  on  deformation  and  plastic 
strain  needs  to  be  evaluated. 

6.2  Flow  Computations 

6.2.1  Flow  Over  Flat  Plate 

To  validate  the  viscous  flow  code  and  to  assess  mathematical  models  for  artificial  dissipation 
the  problem  of  supersonic  flow  over  a  flat  plate  was  solved.  This  problem  has  been  solved 
by  several  investigators  starting  with  Carter  (NASA  TR  R-385)  who  obtained  a  numerical 
solution  using  a  finite  difference  scheme.  The  problem  statement  is  shown  in  Fig.  19a. 
The  mesh  from  adaptive  refinement  of  the  finite  model  is  shown  superimposed  on  density 
contours  in  Fig.  19b.  Temperature  profiles  at  the  exit  plane  are  shown  in  Fig.  19c. 
The  temperature  profiles  shown  are,  in  general,  not  representative  of  high  speed  flows 
because  of  the  high  wall  temperature  boundary  condition.  In  reviewing  the  results  from 
several  different  artificial  dissipation  models,  the  McCormack-Baldwin  model  gives  slightly 
superior  results  to  those  obtained  with  Lapidus  dissipation,  but  it  is  more  expensive  since 
it  requires  extraction  of  second  derivatives  of  the  pressure.  The  successful  solution  of 
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the  Carter  problem  demonstrates  that  the  basic  algorithm  is  working  effectively,  but  work 
remains  to  be  done  to  validate  the  code  for  higher  Mach  number  flows  with  stronger  viscous 
effects. 


6.2.2  Boundary  Layer  Flow 

To  further  evaluate  the  viscous  flow  code,  a  high  speed  boundary  layer  flow  over  a  flat 
plate  is  currently  being  analyzed.  The  problem  statement  is  presented  in  Fig.  20a.  The 
temperature  profile  at  the  inflow  plant  is  shown  in  Fig.  20b.  The  profile  shown  is  repre¬ 
sentative  of  high  speed  flows  and  must  be  predicted  accurately  for  wall  heating  rates  to  be 
predicted  well.  An  important  role  of  this  problem  is  to  provide  insight  into  the  effect  of 
artificial  dissipation  on  the  temperature  distribution  and  the  wall  heating  rates.  Artificial 
dissipation  tends  to  cause  boundary  layers  to  be  too  thick,  distorts  temperature  profiles 
and  can  adversely  affect  wall  heating  rate  predictions.  The  computed  temperature  profile 
at  the  outflow  plane  in  Fig.  20b  shows  only  small  effects  of  artificial  dissipation ,  but  heat¬ 
ing  rates  (not  shown)  computed  from  this  solution  are  below  the  expected  values.  Work 
on  this  problem  is  in  progress.  The  results  so  far  indicate  that  as  the  level  of  artificial 
dissipation  is  reduced,  the  wall  heating  rates  will  improve  in  accuracy.  Successful  comple¬ 
tion  of  this  problem  will  mean  that  the  integration  of  the  flow  analysis  with  a  deformed 
structure  can  be  initiated. 
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7  Research  Forecast 


In  this  initial  contract  year,  the  basic  computational  schemes  for  solving  a  strong  F-T-S 
interaction  problem  have  been  developed  assuming  structured  temperatures  are  below  the 
melting  point.  In  the  next  year,  a  major  goal  will  be  to  integrate  the  components  to  solve 
a  problem  with  a  complete  F-T-S  interaction.  Specific  objectives  for  the  next  year  are: 

1.  Further  enhancement  of  the  thermo-viscoplasticity  code  to  provide: 

(a)  adaptive  mesh  refinement 

(b)  inclusion  of  state  variables  in  the  constitutive  model  to  incorporate  material 
damage 

2.  Further  development  of  the  viscous  flow  code  to  predict  aerodynamic  heating  rates, 
surface  pressures  and  shearing  stresses. 

3.  Development  of  a  mesh  movement  algorithm  for  interfacing  flow  meshes  with  the 
mesh  for  the  deformed  structure. 

4.  Numerical  solutions  for: 

(a)  thermo-viscoplastic  response  of  a  convectively  cooled  structure  with  damage 
predictions  for  multiple  loading  cycles. 

(b)  more  difficult  viscous  flow  problems  with  flow  recirculations 

5.  Development  of  thermal-structural  formulations  for  metallic  materials  where  tem¬ 
peratures  can  exceed  the  melting  temperature. 

6.  Presentation  of  a  paper  at  the  30th  AIAA  Structures,  Structural  Dynamics,  and 
Materials  Conference  to  be  held  April  3-5, 1989  in  Mobile,  Alabama. 
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EQUILIBRIUM  SURFACE  TEMPERATURES 

M£D=8;  q  =  2200  psf 


T,  °F 
3260 


Estimated  National  Aerospace  Plane  Surface  Temperatures  (NASA  TM  S9143) 


Figure  2:  Fluid-thermal-stmctural  Interactions  on  an  Aerospace  Plane  Scranyet  Engine 
Leading  Edge 
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a)  Aerodynamic  Heating  on  Undeformed  Structure 
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b)  Plastic  Deformation  Induced  by  High  Aerodynamic  Heating  on  De¬ 
formed  Structure 


Figure  3;  FTS  Interactions  in  Hypersonic  Flight 
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Figure  4:  Convectively-caoleci  Structure 
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Uniform  Fluid 
Velocity  Profile 


Constant  Mass  Flow 


Figure  5;  Engineering  Model  of  Coolant  Passage  Heat  Transfer 


Figure  6:  Finite  Element  Model  of  Coolant  Passage 
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Figure  8:  Flow  Chart  for  Adaptive  Ti: 


cooling  (665  B) 


(a)  Thermal-structural  model 


Figure  9:  \D  Model  of  Convectively  Cooled  Structure 


(a)  Thermal-structural  model 


(b)  Convection  coefficient  history 


Figure  11:  2D  Model  of  Convectively  Cooled  Structure 


(b)  Temperature  contours  at  t  =  0.5s 


Figure  12:  Thermal  Response  of  2D  Model 
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Figure  19:  Supersonic  Viscous  Flow  Over  Flat  Plate  (Carter  Problem) 
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(b)  Finite  element  model  and  density  contours 


(c)  Exit  temperature  profiles 


Figure  10:  Supersonic  Viscous  Flow  Over  Fiat  Plate  (Carter  Problem) 
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(a)  Problem  statement 
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outflow  profile 


(b)  Temperature  profiles 


Figure  20:  Supersonic  Boundary 


Table  1 :  Data  for  Thermal  Problem 


1 .  Skin  and  primary  structure 


p  =  0.283  lb* /in3 


Temperature  Dependent  Conductivity  and  Specific  Heat 


T  [°RJ  0  500  1500  3000 

K  [BTU  /  in-s-R]  1.23  x  KH  1.8  xlO*4  3.25  xlO-4  6.0  xlO*4 

Cp  [BTU/HvRJ  0.12  0.135  0.177  0.26 


2.  Coolant 

m  =7.16  10  4  lbm  /  sec 

Cp  =  3.48  BTU  11^-  R 

T0  *  665  R 

h  «6.2xl0-3  BTU/s-in‘-R 

3.  Surface  convection 

T„  *  5345  R 

b  *  6.0xl0‘5  BTu/s-iiir-R 

Additional  b  attimc0£t£0.5sce;  ’ 

OS  x<  0.1  ;  h*u  =  0 

O.lSxSO.2:  *3.0xl(Peos  (JLdUil  a  BTU/s-in2-R 

0.05  2  ' 

02  <  x  £  0.4 :  h*ya0 

N 

•  ’  '  v, 

4.  Surface  radiation 

0  =  3.30633  x  10*  B1U  /  s  -  in2  ~  R4 
£  «  0.8 
a  *0.8 
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Table  2:  Data  for  Viscoplastic  Analysis 


Model  Constants  for  superalloy  B  1900  +  Hf 


1.  Temperature  -  independent  constants 

D0  =  1.0  x  104  sec  4 

Zj  =  3000  MPa 

Z3  =  1150  MPa 

rl  -  2 

r2  =  2 

n^  =  0  MPa4 

m2  =  0 

2.  Temperature  -  Dependent  Parameters 


T 

[°C] 

<760 

n 

1.055 

Z2  =  Z0 

[MPa] 

2700 

Ai 

[sec-1] 

0 

a2 

[sec-1] 

0 

Elastic  constants 

871 

982 

1093 

1.03 

0.850 

0.70 

2400 

1900 

1200 

0.0055 

0.02 

0.25 

0 

0 

0 

E  =  1.987  x  105  +  16.78  T  -  0.1034  T2  +  1.143  x  10-5  T5  fMPaJ 

G  =  8.650  x  104  -  17.58  +  0.02321  T2  -  3.464  x  10*5  T3  [MPa] 


(with  temperature  in  °C) 
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